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Abstract. We derive effective interaction potentials between hard, spherical colloidal 
particles and star-branched polyelectrolytes of various functionalities / and smaller 
size than the colloids. The effective interactions are based on a Derjaguin-likc 
approximation, which is based on previously derived potentials acting between 
polyelectrolyte stars and planar walls. On the basis of these interactions we 
subsequently calculate the demixing binodals of the binary colloid-polyelectrolyte star 
mixture, employing standard tools from liquid-state theory. We find that the mixture 
is indeed unstable at moderately high overall concentrations. The system becomes 
more unstable with respect to demixing as the star functionality and the size ratio 
grow. 

PACS numbers: 61.20.-p, 61.20.Gy, 64.70.-p 
1. Introduction 

Polyelectrolyte stars (PE's) are complex macromolecules that have attracted a lot of 
interest in the recent past. They consist of / polymer chains, all attached on a common 
centre, and carrying ionizable groups along their backbones. Solution of these molecules 
in a polar solvent results into dissociation of the groups, so that the chains turn into 
polyelectrolytes and stretch considerably with respect to their neutral conformations. 
Already in the early 1990s, the importance of the stretched PE-chains in stabilizing 
colloidal suspensions has been pointed out and analysed theoretically by Pincus [1] 
employing scaling theory as well as, more recently, by Wang and Denton [2] using 
linear-response theory. A distinguishing feature of PE-stars is their ability to adsorb 
the vast majority of the released counterions into their interior, creating thereby an 
inhomogeneous cloud of entropically trapped particles that provides a strong entropic 
barrier against coagulation [H, [2], [31 SI El El d El E] - The development of accurate effective 
interactions between the PE-stars 0, 13 El E] has led to predictions regarding their overall 
phase behaviour with emphasis on crystallization [101 [11], which has recently received 
experimental corroboration [12] . 
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Though a great deal has thus been learned regarding the behaviour of one- 
component solutions of PE-stars, the question of the influence of these ultrasoft colloids 
on solutions of hard colloids has not been investigated thus far. At the same time, 
the behaviour of PE-stars in the vicinity of planar or curved hard surfaces (such as 
a larger colloidal particle) is an issue of considerable interest, due to the possibility 
of manipulating the conformation of the PE-star by suitably changing the surface's 
geometry or physical characteristics |T3l O EE] ■ Recently, the properties of PE-stars 
close to hard, planar walls were investigated in detail by means of computer simulations 
and theory [16] . It has been found that the geometrical constraint of the planar wall does 
not affect the ability of the PE-stars to absorb the vast majority of their counterions. 
In addition, a new mechanism giving rise to a wall-star repulsion has been discovered, 
which rests on compression of stiff star chains against the neighboring wall. In this 
work, we proceed to the full, many body problem of a collection of PE-stars and neutral 
colloids, which can be seen as curved walls. Basing on the results of Ref. [16J, we 
investigate the structure of the mixture and find that it is unstable against demixing as 
the concentration becomes sufficiently high. This work serves, thereby, as the reference 
point for future investigations on the effects of adding charge to the colloidal particles. 
It is specular to recently published work on mixtures of charged colloids with uncharged 
polymers [T7j, since in our case the colloids are neutral and the (star-branched) polymers 
are charged. 

The rest of this paper is organised as follows: in sec. [2] we introduce the colloid- 
colloid and PE-star-PE-star effective interactions and we derive the cross interaction, 
based on previous results on the PE-star interaction potential with a planar wall. In sec. 
[3]we present our method for calculating structure and thermodynamics by employing the 
aforementioned interactions in combination with two-component liquid integral equation 
theories. In sec. H]we present our results for various regimes of the parameter space as 
well as the overall phase diagrams of the mixture. Finally, in sec. [5] we summarize and 
draw our conclusions. 

2. Effective pair potentials 

The system under investigation is a binary colloid-PE-star mixture. The colloids 
are coded with the subscript 'c' and the PE-stars with 's'. The mixture contains, 
thus, N c spherical, neutral colloids with diameter a c (radius R c ) and iV s PE-stars in 
aqueous solution. The stars can be characterised by their degree of polymerization N p , 
functionality /, and charging fraction a. Thereby, the / chains of each star are charged 
in a periodical manner in such a way that every (1 / a)-th monomer carries a charge e. As 
a result, every star carries a total bare charge Q = eafN p , leaving behind M = afN p 
monovalent, oppositely charged counterions in the mixture due to the requirement that 
the system must remain electro-neutral as a whole. With er s referring to the stars' 
diameter, i.e., twice the average centre-to-end distance R s of the arms, we define the 
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size ratio q between the two species as 



Q = o- s /a c . 



(1) 



Within this work, we will only consider PE-stars that are smaller than the colloids, hence 
q < 1. The degree of polymerization of every arm, N p , and the charging ratio a play 
a crucial role because they determine the number of released counterions M mentioned 
above. The latter are, in turn, mainly responsible for the emergence of the star-star 
[H El [8] and the star-colloid effective repulsions [16j, due to the loss of entropy they 
experience when two such objects approach close to each other, see also eq. ([3]) in what 
follows. In this work, we fix N p = 50 and a = 1/3 throughout. Generalizations to other 
values of a and N p can follow by appropriately taking into account the dependence of 
M on these parameters. Thereby, the two remaining single-molecule parameters that 
we vary are the stars' functionality / and the size ratio q. 

The thermodynamic parameters are the partial number densities pi = Ni/V 
(i = c, s) of the respective species and the absolute temperature T. Alternatively, we can 
work with the concentrations X{ = Ni/N and the total number density p = N/V, with 
the total particle number N = N c + N s in the overall volume V of our model system. We 
will consider constant, room temperature (T = 300 K) throughout this work. This is the 
temperature for which the star-star effective interactions [7J [8] and the PE-star-planar 
wall potentials [16] have been derived, based on the value Ab = 7.1 A for the Bjerrum 
length in aqueous solvents. As usual, we define the inverse thermal energy (3 = 1/(/cbT), 
with &b denoting Boltzmann's constant. 

The starting point for all considerations to follow are the effective pair potentials 
between the constituent mesoscopic particles, having integrated out all the monomer, 
solvent and counterions degrees of freedom. When introducing this set of interactions as 
an input quantity into the full two-component integral equation theory described in more 
detail in Sec. [3l we can in principle completely access the structure and thermodynamics 
of the system at hand. 

2.1. The colloid-colloid and PE-star-PE-star interactions 

The effective colloid-colloid interaction at centre-to-centre distance r is simply taken to 
be a pure hard sphere (HS) potential, namely: 



A lot of work concerning effective PE-star-PE-star interactions was done in 
the recent past by Jusufi and co-workers [7J [S]. They employed monomer- resolved 
Molecular Dynamics (MD) simulations and analytical theories and found an ultra-soft, 
bounded, density-dependent effective interaction governed by the entropic repulsions of 
counterions trapped in the interior of the stars. The good agreement between simulations 
and theory even allowed them to put forward analytic expressions for the full pair 
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potential at arbitrary star separations. The effective potential has a weak density- 
dependence, which however disappears when the star density exceeds its overlap value 
p*. In this case, all counterions are absorbed within the stars, whose bare charges are 
therefore completely compensated. Thus, the effective potential vanishes identically 
for centre-to-centre distances r > a s . For overlapping distances r < cr s , there is no 
dependence on the concentration anymore and only the trapped counterions' entropy 
contributes to the star-star interaction, for this reason reading for r < o~ s = qo c as 

(3V ss {r) , I N 2 I ^ ln2 fe) , (N 2 
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In eq. (j3J) above, N 2 is the number of spherically trapped counterions of a single star. It 
does not coincide with the number of released counterions, M, because the number N\ 
of Manning-condensed counterions [19] does not contribute to the effective interaction 
and must be excluded: thus N 2 = M — N\. Extensive simulations [H [H [161 EI] have 
shown that the relative population of counterions in the two possible states is essentially 
independent of r. Thus, we fix Ni to the value measured in MD simulations made during 
the investigation of PE-stars in planar confinement [16j. The fraction N±/M typically 
grows with increasing a and covers ranges between 30% and 50%. 

Clearly, the interaction V ss (r) of eq. ([3]) vanishes, along with its first derivative with 
respect to r, at r = qo~ c , guaranteeing the smooth transition to the region r > qcr c , in 
which V ss (r) = 0. The latter feature is, strictly speaking, valid only for star densities 
exceeding the overlap value p* [TT]. For p s < p*, a Yukawa tail exists, emerging from 
the Coulomb interaction between the non-neutralised PE-stars and screened by the free 
counterions. For the purposes of simplicity, we ignore this small contribution, because 
the number of released counterions from multiarm PE-stars constitutes, at all densities, 
a tiny fraction of the total number of counterions M [8] , as confirmed by the very small 
values of experimentally measured osmotic coefficients of PE-star solutions [IE] • 



2.2. The cross interaction 

In order to complete the set of effective pair potentials needed to describe the binary 
mixture within the framework of a full two-component picture, we have to specify the 
colloid-PE-star cross interaction. Thereby, we proceed along the lines of Ref. [22] to 
derive the desired potential for small q- values based on results for the effective repulsion 
in the case where a PE-star is brought within a distance z from a hard, flat wall [16, 21J. 

To begin with, let V sw (z) be the star-wall interaction and F sw (z) = —dV sw (z)/dz 
the corresponding force for a PE-star with all its counterions absorbed, i.e., for densities 
p s beyond the overlap density (see also previous Sec. I2.ip . Then, for the geometry shown 
in fig. [U^a), the force is related to the osmotic pressure Il(s) exerted by the star on the 
surface of the wall via integration of the normal component of the latter along the area 




of contact pQ: 

/OO /*OG 
dy yU(s) cos ^ = 2nz J dsU(s). (4) 

Using the above eq. HI we can directly obtain the functional form for the osmotic pressure 
II(z), provided that the functional form for the star-wall force F sw (z) is known: 



2vrd2 

The same ideas can in principle be applied for a PE-star in the vicinity of a spherical, 
hard colloid, i.e., a hard sphere. Again, integrating the osmotic pressure along the area 
of contact between both objects yields the force acting on the centres of the mesoscopic 
particles. Pursuant to the geometry of the problem, see fig. DJb), and paying regard to 
the underlying symmetry, we get as result for the colloid-PE-star cross force F* s (z): 

F* s ( z ) = —^ d9 sin#n(s)costf. (6) 

2 Jo 

Here, the upper integration boundary 6* max can be acquired by the condition that IT(s) 
must vanish identically for all 9 > 9 mSuX . It is possible to eliminate the polar angles i? 
and 9 emanating from the centres of the PE-star and the colloid, respectively, in favor 
of the distance s between the star centre and the point on the colloid's surface that 
is determined by the aforementioned angles. In doing so, we use geometrical relations 
evident from the sketch in fig. QJb), and finally obtain: 

Tin f s xnax 

FXz) = - c 2 / ds [(a c + 2zf - al + As 2 } U(s). (7) 
I (a c + 2z) Jz 

Again, we may obtain the maximum integration distance s max (without any need to 
calculate 9 milx before) simply by demanding that II (s) must be equal to zero for all 
s > s max . For such values of s, the integrand as a whole obviously vanishes and there 
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are no contributions to the result of the integration anymore. Presumed the functional 
form for the osmotic pressure is known, such identification of s max is rather easily feasible. 

Since we want to consider small values q < 0.3 of the size ratio only, the stars discern 
the colloidal surface to be rather weakly bent compared to a flat wall, i.e., the radius 
of curvature is large in terms of the star diameter a s . Therefore, it is a reasonable 
approximation to assume that the osmotic pressure remains almost unchanged with 
respect to the situation where a PE-star is brought in contact with a planar wall. 
Consequently, we may combine eqs. ([5]) and (JTj) to obtain a sound estimate for the 
effective force F* s (z) as a function of distance of the star centre and the colloid's surface. 
Note that in our special case s max is of the order of the star radius i? s . This fact becomes 
evident from eq. ([5]) if one takes into account that the typical range for the star-wall 
force F sw (z) is also approximately i? s or at the utmost slightly bigger due to effects of a 
chain compression at the hard wall (see below) [16]. Clearly, the corresponding potential 
is received by a simple, one-dimensional integration: 



In fig. [2] we show the shape of V* s (z) for q = 0.2 and different values of the stars' 
functionality /. In order to demonstrate the importance of the so-called compression 
term adding to the star-wall force F sw besides electrostatic-entropic contributions 
[T6] . we additionally included colloid-star potentials calculated on the basis of the 
electrostatic and entropic star-wall forces alone. Since there are striking deviations, we 
can clearly expect such devolved compression effects to influence the phase behaviour 
of the mixture. 

Finally, we need to express the effective potential as a function of the particles' 
centre-to-centre separation r instead of the centre-to-surface distance z. Thereby, we 
have to take into account that the star centre is strictly forbidden to penetrate the 
volume of the colloid. Thus, the total cross interaction features a hard core plus the 
soft, purely repulsive tail as obtained from the above eq. ([8]) and finally writes as: 



3. Determination of the structure and thermodynamics of the mixture 

In this section, we describe the basic principles of liquid integral equation theory for 
binary mixtures^ and how to subsequently access the thermodynamics of the system. 
In general, the pair structure of the system at hand (and analogously any other two- 
component system) is fully described by three independent total correlation functions 
hij(r) with i,j = c, s. Hereby, we already allowed for the symmetry with respect to 

§ A further generalization of the theoretical approach from v = 2 to v > 2 components in the mixture 
is straightforward. But since we are only interested in binary systems within the framework of this 
paper, we limit ourselves to that special case in order to keep the delineation as concise as possible. 




(8) 




(9) 
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Figure 2. Effective colloid-PE-star potentials with and without provision for 
compression effects as a function of the centre-to-surface distance z. Here, we have 
chosen the parameters q — 0.2 and (a) / = 10, (b) / = 18, and (c) / = 30. In the 
legend boxes, the numbers of condensed counterions N\ used as fit parameters in Ref. 
[I6j are specified for sake of completeness. 



exchange of the indices, i.e., hij(r) = hji(r). Closely related to the total correlation 
functions are the so-called direct correlation functions (def's) c^(r). Following the same 
symmetry argument again, there exist only three independent def 's. In what follows, we 
will denote the Fourier transforms of hij(r) and Cy(r) as hij(r) and Cjj(r), respectively. 

The above-mentioned connection between the total and direct correlation functions 
is quantitatively incorporated via the multicomponent generalization of the well- 
known and commonly used Ornstein-Zernike (OZ) relation, which in its Fourier space 
representation reads as [231 I2U 125] : 

H(Jfe) = C(ife) + C(Jfe)-D-H(Jfe). (10) 

Here, H(fc) and C(k) are symmetric (2 x 2) matrices whose elements are constituted 
by the total and direct correlation functions, respectively, and D is a diagonal (2 x 2) 
matrix containing the partial densities characterising the composition of the system 
under investigation, i.e., 

H(fc)l =kj(k), (11) 
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C(A;) =Cij(k), (12) 

L J ij 

IP]« = M*- (13) 
Evidently, eq. (fTUl) can be rewritten yielding the equivalent matrix relation 

H(fc) = [l -C(fc)-Dl _1 -C(jfc), (14) 

r i _1 

with the identity matrix 1 and the matrix inverse 1 — C(fc) • D . Defining A(k) = 
PsPc[css(k)c cc (k) — c^(k)] and E{k) = p s c ss (k) + p c c cc {k) and returning to a component- 
by-component notation, the latter can consistently be expressed in the following fashion: 

h (k) ^ ~ ' A(fe) ' 5ij (15) 
KAk) ~ 1 + A(k)-E(k) ' [ } 

The linear algebraic system of eq. (IT51) . provides three independent equations 
coupling six yet unknown functions hij(k) and Cij(k). In order to completely determine 
that set of functions, we therefore need to supply three additional relations to close 
and subsequently solve the system of equations. There are several popular choices 
for these so-called closures, e.g., the Percus-Yevick (PY) or hypernetted-chain (HNC) 
approximations in their respective two-component generalizations. While the former is 
known to generate reliable results for short-range interactions, the latter furnishes very 
accurate estimates for the pair structure in case of long-ranged, soft potentials. Neither 
the PY nor the HNC closure are thermodynamically consistent, however, and in our case 
this is a crucial factor, since we are interested in the calculation of phase boundaries, 
which should not depend on the route chosen to calculate the free energies. Thus, we 
resort to the Rogers- Young (RY) closure [26] , in which thermodynamic consistency can 
be enforced. In its multicomponent version the RY-closure reads as: 

gij {r) = exp [-pVy] ■ 1 1 + J -j-^ 1 , (16) 

where gij(r) = hij(r)+l are the so-called radial distribution functions and we introduced 
new auxiliary functions Xij{ r ) = hij( r ) ~ c ij{ r )- Vij{ r ) refers to the pair interactions 
between species i and j as presented in sec. |2j It may be again emphasised that the main 
benefit we gain from using the modified relation (jT6l) is closely related to the hybrid 
character of the latter. Due to the fact that any closure constitutes an approximation, we 
in general obtain different results for the partial and total isothermal compressibilities as 
calculated via either the virial or the fluctuation route (see below), as already mentioned 
above. But the three mixing functions emerging in eq. ( |T6l) above and given by 

fij( r ) = 1 - ex P [-Cijr] , (17) 

with being the so-called self-consistency parameters, now allow us to address 
this problem and to appropriately match the isothermal compressibilities. Since it 
is sufficient to apply a single consistency condition only, namely the requirement of 
equality of the system's total virial and fluctuation isothermal compressibilities, the 
usual approach is to employ just one individual parameter Qj = ( for all components. 
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Hence, only a single mixing function fij(r) = f(r) remains. However, multi-parameter 
versions of the RY closure have nevertheless also been proposed some years ago [27J, 
accordingly demanding the equality of all the partial compressibilities. It is easy to check 
that for C = and ( = oo the multicomponent PY and HNC closures, respectively, are 
recovered from eq. (TIBjfjj l. 

Now, we have to address in more detail the issue of calculating the total isothermal 
compressibility following the different routes. At first, we concern ourselves with the 
virial compressibility Kj.. The total pressure P of the system at hand, including both 
ideal and excess contributions, takes the form |25j : 

2,71 0^ f°° 

PP = p--£-Y,Y l x i x j / drr 3 V;Ar) giJ (r), (18) 



with V(j{r) = —dVij(r)/dr being the different pair potentials' derivatives with respect 
to the inter-particle distance r. Provided the pressure pursuant to eq. (ITS]) is known, 
can be obtained by differentiating with respect to the total density p while the partial 
concentrations Xi £1X6 kept fixed: 

-i 



pk^Tn] 



~d{(3P) 




dp 





(19) 



In order to evaluate the fluctuation compressibility Kj,, we initially introduce the 
three partial structure factors Sij(k). As the correlation functions and the radial 
distribution functions, respectively, they also describe the structure of the system: 

Sij(k) = Sij + ^/~p~p~hij{k). (20) 

While for the one-component case the compressibility can simply be obtained as the 
(k = 0)-value of the static structure factor, i.e., S(k = 0) = pksTK^, things are 
a bit more complicated for binary (or multicomponent, u > 2) mixtures. Here, in 
generalization of the one-component situation, the compressibility can finally be written 
using the following expression [28, 29, 30J: 



, T a _ S ss (0)S cc (0)-ff c 2 s (0) (m 

9 B T x c S ss (0) + x s S cc (0) - 2^x-x-Sl(0) ■ { ) 

Based on the knowledge of the partial correlation functions hij(r) and structure 
factors Sij(k) as obtained by (numerically) solving the OZ relation, eq. (JTU1) . and using 
the RY closure, eq. f|T6|) . we can in principle completely access the thermodynamics of 
the system at hand. In order to calculate the binodal lines, a very convenient quantity 
to consider is the concentration structure factor S con (k). It is a linear combination of 

| When using the RY closure, the correlation functions obviously, besides their inherent density 
dependence, parametrically depend on the mixing parameter £, i.e., /ly = hij(r; p c , p s ,C) an d 
Cij — Cij (r; p c , p s , £)■ The same must obviously hold for all quantities deduced from these two functions. 
Note that we will nevertheless throughout this paper drop both the pj's and £ from the respective 
parameter lists in order not to overcrowd our notation. 
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all the partial structure factors, whereby the corresponding pre-factors are determined 
by the different species' concentrations Xj, namely: 

*S'con(^) = x c x s S cc (k) + x s x c S ss (k) — 2(x c x s ) 3 / S cs (k). (22) 

Now, let P be the total pressure according to the above eq. ({IB]) and g(x s ,P,T) = 
G(x s , N, P,T)/N the Gibbs free energy G(x B , N, P,T) per particle. Then, the second 
derivative of the former is connected to the concentration structure factor S con (k) by 
means of the sum rule [31j [321 [33] : 

where we have added the concentration x s as a second argument to S con (k) to emphasise 
this dependence. This differential equation has to be integrated along an isobar for any 
prescribed value of the pressure P* = PPcr^ = const, to obtain the Gibbs free energy 
from the structural data, S con (k = 0;x s ). A detailed analysis of the limiting behaviour 
of g"(x s ) shows a divergence as l/x s for x s — > and as 1/(1 — x s ) for x s — > 1 |33j. In 
order to avoid any technical difficulties when numerically integrating, we a priori split 
the Gibbs free energy g(x s ) into a term that arises from its ideal part and a remainder, 
which we call excess partial, g cx (x s ): 

Pg(x s ) = (1 - x s ) hi(l - x s ) + x s ln(x s ) 

+ 3(1 - x s ) ln(A c /a c ) + 3x s ln(A s /a c ) + Pg^x,), (24) 

with the thermal de Broglie wavelengths A c s of the colloids and the stars, respectively. 
Taking the second derivative in the above eq. (I24|) again, we obtain: 

/3gT(x a ) = - + T ^— + C^). (25) 

x s i x s 

Thus, the ideal part of the Gibbs free energy is exclusively responsible for the appearance 
of the aforementioned divergences at the integration boundaries and the modified 
second-order differential equation 

*S'con(0, X s ) 1 "^8/ 

for the excess Gibbs free energy alone is obviously free of any diverging terms. We can 
therefore easily solve it numerically. Subsequent addition of the analytically known ideal 
term gid(^s) = (1 — x s ) ln(l — x s ) + a; s ln(a; s ) directly yields the total Gibbs free energy 
per particle that we are interested in. The two terms involving the thermal de Broglie 
wavelength are linear in x s ; they only provide a shifting of the chemical potentials and 
can be dropped. 

Thermodynamic stability requires that g(x s ) is convex [34]. In case we encounter 
some x s -region where g"(x s ) < the binary mixture features a fluid-fluid demixing 
transition. In that sense, eqs. (|23|) and (|26l) . respectively, allow us to investigate the 

The 'excess' part g e x(x s ) in eq. (|24[) . includes a term ln(pcr^) that arises from the original ideal part 
but which does not cause any divergences at the limits x s — > and x s — ► 1, which we seek to remove. 
Thus we readsorb it into a redefined excess part, which can be integrated without problems. 
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thermodynamics and the phase behaviour of the system at hand by providing a tool to 
compute the Gibbs free energy (per particle). The corresponding phase boundaries can 
be calculated using Maxwell's common tangent construction, which guarantees that the 
chemical potentials, the same between both coexisting phases. Since we are in 

a situation where we moreover fixed the pressure P* of the mixture and its absolute 
temperature T, all conditions for phase coexistence are clearly fulfilled. Concretely, the 
common tangent construction amounts to solving the coupled equations 

9W = 9'(xf) (27) 

and 

9ti) ~ xW(xl) = 9(xf) - xfg'(xf) (28) 

for the concentrations xl ,n of the coexisting phases I and II. 

In integrating eq. (1251) above and adding the ideal terms, one obtains the Gibbs free 
energy per particle, g(x s ) modulo an undetermined linear function C\x s + Co with the 
constants G\ and C to be fixed by appropriate boundary conditions. As is clear from 
eqs. (127]) and (|28|) above, such a linear term is anyway immaterial from the determination 
of phase boundaries and, in practice, it can be ignored on the same grounds that the 
terms involving the thermal de Broglie wavelengths in eq. (T2~^]) have been dropped. 
Nevertheless, the constants C\ and Co can be determined as follows. Taking into account 
that the Gibbs free energy G(N, P,x s ,T) is an extensive function but in its natural 
argument list there is only one extensive variable, namely the number of particles N, 
Euler's theorem [35] asserts the function g to have the form: 

g(x s ) = x s p s (x s ) + (1 - x s )p c (x s ). (29) 

For both limiting one-component cases, i.e., if no stars (x 8 = 0) or no colloids (x s = 1) 
are present in the system, the following relation holds true: 

Pg = f + pf = Hp) + Ax + pfL (so) 

where / = (3F/N denotes the Helmholtz free energy per particle and /' its derivative 
with respect to the density p\ the subscript 'ex' refers to the excess part of /. On the 
other hand, f' ex is connected to the excess pressure P ex [as known from eq. ( fl8|) above] 
via the equation 

f' ex = (3P ex /p 2 . (31) 

Following eq. ( 13T1) . f^ x can be obtained by integrating the ratio P cx / p 2 with respect 
to p and applying the additional boundary condition f ex (p — > 0) = 0. Once the 
Helmholtz free energies for the pure colloid and PE-star systems are known this way, the 
corresponding chemical potentials /x c (0) and /i s (l) can be calculated and the conditions 
g(0) = p c (0) and g(l) = p s (l) for any arbitrary pressure P [cf. eqs. (|29l) and (1301) 
above], yield Co and C\. Note that for the pure colloidal system, x s = 0, we can avoid 
the integration route to compute the pressure, by using the accurate Carnahan-Starling 
expressions for hard-spheres [36], which also turn out to be consistent with the one 
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Figure 3. Examples of the second derivative of the Gibbs free energy per particle, 
g"(x s ), plotted against the star concentration x s for stars with functionality / = 30, 
PE-star-colloid size ratio q = 0.2, and two different pressures P* = /3Ptj%- Symbols are 
calculated from the OZ relation, lines were obtained by cubic spline interpolation. Note 
that the a; s -interval where we are not able to numerically solve the integral equations 
grows distinctly upon increasing the pressure, i.e., as we move away from the critical 
point. Moreover, the plot illustrates the limiting behaviour of g"(x s ) as l/x s for x s — * 
and as 1/(1 — x s ) for the opposite case x s — ► 1, respectively. 

calculated from the Rogers- Young route, based on eq. (JTBj) and our results for the radial 
distribution function g(r). 

Note that, when crossing the spinodal line in the density plane, the long wavelength 
limits of the partial structure factors, S{j(k — > 0), take non-physical values. This 
behaviour expresses the system's physical instability against a possible fluid-fluid phase 
separation. Thereto, it is not feasible to (numerically) solve the integral equations 
anymore once we reached the spinodal; in fact, integral equation theories themselves 
break down before the spinodal is reached, yet after the binodal [37] . Consequently, 
depending on the total pressure and above a certain threshold value of the same, 
P > -Pthr, the concentration structure factor Scon(0',x s ) is unknown over some interval 
Ax s (P). Hence, we need to appropriately interpolate S con (0; x s ) in order to obtain the 
second derivative of the Gibbs free energy per particle for all < x s < 1 and, in this 
way, to allow for the integration of the differential equations ( 123]) or (I26I) . respectively. 
Along the lines of Ref. |33j, we perform this necessary interpolation using cubic splines. 
In order to illustrate the whole procedure, fig. [3] shows the function g"(x s ) as computed 
from the OZ equation together with its cubic spline interpolation for a representative 
parameter combination and two different pressures. Moreover, in fig. HI we plotted the 
corresponding Gibbs free energy g(x s ) for the lower one of these pressures. In addition, 
the inset depicts Maxwell's common tangent construction used to compute the star 





Figure 4. Gibbs free energy per particle g(x s ) vs. the star concentration x s , plotted 
for / = 30, q = 0.2, and P* = 50.0. The curve was obtained via integrating the 
interpolated function g"(x s ) twice according to the procedure delineated in the main 
text, whereby we subtracted an arbitrary linear function afterwards. The inset shows 
g(x s ) with an differently scaled £ s -axis in order to highlight the concave parts of the 
function and to show Maxwell's common tangent construction. 



concentrations for the coexisting phases. 

It may be emphasised that the results for the binodal do not depend on the concrete 
interpolation scheme, at least as long as the numerical methods used to solve the OZ 
relation are able to precisely reach the spinodal, i.e., the points where the structure 
factors diverge for k — *■ 0. Admittedly, this is not always strictly the case since the 
numerical schemes we employed to calculate correlation functions and corresponding 
structure factors, respectively, may break down slightly before the spinodal is reached. 
Accordingly, small inaccuracies induced by the interpolation procedure arise which grow 
with increasing width of the gap region Ax s (P), or to put it in other words, with 
increasing pressure P, i.e., if we move away from the critical point. As long as the 
aforementioned interval where no solution of the integral equations can be found is 
rather small, we expect the interpolation to be reliable, while for higher pressures the 
received binodals are of more approximate character. Nevertheless, they still show a 
very reasonable behaviour. We are going to discuss the results for the phase diagrams 
in detail in sec. 14.31 
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4. Results 

4-1. Low colloid-density limit 

Based on the radial distribution functions g(r) as obtained by the OZ relation closed 
with the RY closure, eqs. ( fTOi) and (TTBI) . we may map our two-component mixture onto 
an effective one-component system of the colloids alone. In doing so, the PE-stars are 
completely traced out, resulting into an effective colloid-colloid interaction where the 
pure hard-sphere potential is masked by additional depletion contributions originating 
in the presence of the stars and the forces they exert on the colloids. To put it in 
other words, the colloid-PE-star interactions cause spatial correlations of the PE-star 
distribution in the vicinity of the colloids, and it is exactly these correlations that 
determine the resulting shape of the depletion potential. Note that the latter in general 
parametrically depends rather on the PE-stars' chemical potential p s or, equivalently, 
the density p r s of a reservoir of stars at the same chemical potential pj g = p s , than on 
their density p s in the real system. Hence, it is in principle more convenient to switch to 
a reservoir representation (p c , p T s ) of the partial densities instead of the original system 
representation (p c , p s ) when considering such effective interactions. Clearly, if the colloid 
density p c takes finite values, it must hold p s 7^ p\. But since we will consider the limiting 
case of low colloid densities p c — > only in what follows, we have p r s = p s again, i.e., 
reservoir and system representation of the partial densities coincide. 

Concretely, the desired mapping^ can be achieved by a so-called inversion of the 
full, two-component results of the integral equations in the low colloid-density limit 
p c — > [331 EH EHl HOI HI]. It can be shown from diagrammatic expansions in the 
framework of the theory of liquids [23J that in this limit the pair correlation function for 
any fluid reduces to the Boltzmann factor g{r) = exp[—f3v{r)]. Here, v{r) denotes the 
pair potential the fluid's constituent particles interact by. According to this relation, the 
effective colloid-colloid potential V e s(r), depending parametrically on both the partial 
colloid and star densities p c and p r s = p s , is obtained as follows: 



Fig. [5] shows examples for the effective colloid-colloid interaction V e s(r) for different 
functionalities / of the stars, partial star densities p T s = p s , and PE-star-colloid size ratios 
q. As one can see from the plots, for distances r > a c the resulting depletion interaction 
mediated by the stars is attractive and features a slightly oscillating behaviour, while 
for inter-particle separations r < a c the bare hard-sphere repulsion remains. In 
particular, fig. [3(a) illustrates that the addition of PE-stars to the mixture results in 

+ Note that the most accurate way to compute effective interactions between two colioidal particles 
in the presence of (smaller) PE-stars is to employ direct computer simulations [33l EH [42], [43] 144] , 
Another way to the depletion potential would in principle be offered by Attard's so-called superposition 
approximation (SA) [45j . But since we want to perform the mapping onto an effective one-component 
system in order to gain some qualitative understanding of the physics of our system only but stick to 
the full two-component picture to quantitatively calculate the binodals of the mixture, we turn down 
such alternative methods within the scope of the paper at hand. 




(32) 
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Figure 5. Effective colloid-colloid depletion potentials V c g(r) as obtained by an 
inversion of the OZ relation. For details concerning the procedure, see main text. We 
have investigated the influence of (a) the partial star density p' s = p s , (b) the PE-star- 
colloid size ratio q, and (c) the stars' functionality / on the functional form of the 
interaction potential. It is evident from the plots that the presence of the stars induces 
an attraction between the colloids in addition to their bare hard-sphere repulsion which 
takes over for distances r < a c . 



both a significant increase of the depth of the attractive potential well and a further 
enhancement of the aforementioned oscillations but does in no way affect the range of 
the attraction. As can be read off from fig. Efb), the latter is determined by the size 
ratio q alone and grows linearly with the diameter of the stars. Furthermore, there is 
a measurable, indeed weak, dependence of the interaction strength on the functionality 
of the PE-stars: The higher the arm number / gets the stronger becomes the effective 
attraction between two colloids, cf. fig. Etc). All these trends are in perfect agreement 
with the common understanding of the physical mechanisms leading to the appearance 
of such an effective attraction: Due to a depletion of the PE-stars in the spatial region 
between a pair of colloids and dependent on the colloids' mutual distance, they are 
hit asymmetrically by the stars from the inside and the outside. Consequently, the 
unbalanced osmotic pressure pushes the colloids together. Clearly, the absolute value 
of this force must grow when increasing the star density p s , simply because there are 
more collisions between PE-stars and colloids. For higher functionalities /, the colloid- 
PE-star cross interaction becomes more repulsive (see sec. 12.21) . i.e., the stars push the 
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colloids harder, thus also leading to a strengthened effective colloid-colloid attraction. 
And finally, the PE-stars' diameter determines whether or not they fit into the spatial 
region between a pair of colloids for a given distance of the two. Hence, the size ratio q 
controls if the stars are expelled from the said region of space, or to put it in other words, 
for what scope of inter-colloidal separations depletion actually takes place. Accordingly, 
the range of the effective force can be altered by changing q. 

The occurrence of oscillations of the effective potential V e g obviously means that 
the attractive minimum is followed by a repulsive barrier whose height is set by the 
concentration of PE-stars in the mixture, see above. In particular, it grows upon 
addition of stars to the system and such behaviour could in case of distinctly high 
and broad maxima in principle lead to micro-phase separation, i.e., cluster formation 
[471 l48l l49l l50| I5T] . But for the physical system we examine and the range of parameters 
we investigate, the barrier remains anyway rather low and narrow. Micro-phase 
separation is therefore not likely to happen. Instead, the type of effective colloid-colloid 
attractions at hand, i.e., an attractive potential valley together with a nearly vanishing 
or at least less-pronounced repulsive barrier, forces the system to develop long-range 
fluctuations upon an increase of the PE-star concentration, consequently favoring the 
possibility of a fluid-fluid demixing transition of the two-component mixture. Such 
behaviour is frequently observed in, e.g., colloid-polymer mixtures [52j [531 [51] . Thus, 
when considering the phase behaviour of our system by calculating its binodals, we 
expect to find evidence for macro-phase separation. This supposition will be endorsed 
by the results of the following section, too. 

4-2. Structure of the mixture 

Before switching over to a presentation of the demixing binodals as obtained via the 
procedure described in detail in Sec. [3] of this paper, i.e., initially calculating the Gibbs 
free energy g(x s ) with both the temperature T and the pressure P kept fixed and 
subsequently identifying the sought-for coexisting fluid phases using Maxwell's common 
tangent construction for the concave parts of that function (see, in particular, figs. [3] 
andH]), it is useful to study partial pair correlation functions gij{r) and corresponding 
structure factors Sij(k) (i,j = c, s) first. Since these quantities completely describe 
the pair structure of the system, we are able to gain detailed insight into the physics 
and phase behaviour of the mixture and to discover, in addition to the findings of the 
previous Sec. 14.14 more evidence that it is reasonable to expect an mixing-demixing 
transition. 

Fig. [6] shows the partial radial distribution functions gij(r) for typical parameters, 
namely a colloid-PE-star mixture with a size ratio of q = 0.3 and the PE-stars having 
/ = 30 arms each. We show results for different mixture compositions, i.e., varying 
partial densities for both species as indicated in the plots. Figs. [6](a) and (b), on 
the one hand, depict the decisive length scales of the problem or, equivalently, the 
typical ranges of the underlying pair potentials as set by the the sizes of colloids 
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Figure 6. Partial radial distribution functions gij(r) (i — c, s) for PE-stars with 
/ = 30 arms, star-colloid size ratio q = 0.3, fixed PE-star density p s a^ — 0.27, and 
colloidal densities of (a) p c a^ = 0.05 or (b) p c a^ — 0.29. For the same value of the 
stars' partial density, the remaining two parts of the figure illustrate the detailed shape 
and the p c -dependence of (c) the star-star correlation function g ss {r) and (d) the cross- 
correlation function g cs (r). For an in-depth discussion and interpretation of the results, 
we refer the reader to the main text. 



and PE-stars, respectively. The distinct height of the colloid-colloid contact value 
9cc{o~ c ) and its further rise upon increasing the PE-star density (not shown in our 
figures) is an obvious manifestation of the mainly attractive character of the effective 
colloid-colloid interactions. In this respect, we again refer the reader to Sec. 14.11 and, 
in particular, eq. (1321) mathematically describing the inversion procedure for the OZ 
relation. On the other hand, when taking a closer look to the whole set of pair 
correlation functions, we find various signs pointing towards the supposable occurrence 
of a demixing transition. The main peaks of both g ss (r) and g cc {r) gain in height when 
adding colloids to the system, while the peak height of the cross- correlation function 
g cs (r) remains essentially the same, see figs. E](a), (b) and (d). In addition, figs. Etc) and 
(d) show an enhancement in the star-star pair correlations and an concurrent depletion 
in the colloid-star correlations for raising colloid densities. The intervals of distances 
affected are remarkably broad, both the range of the enhancement and the depletion 
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Figure 7. Examples of the partial structure factors (a) S cc (k), (b) S sa (k), and (c) 
S cs (k) for PE-star functionality / = 18, size ratio q = 0.2, fixed density of the stars 
PsCTg = 0.12, and several values of the colloidal density /9 c Cc, i.e., different mixture 
compositions. Please note that the line styles in the main plot of part (b) refer to the 
same parameters as explained in the legends of parts (a) and (c) , respectively. The inset 
in (b) addresses a comparison between the colloid-colloid and the star-star structure 
factors for the aforementioned star density and a typical value of the colloid density 
(indicated in the plot) and thereby illustrates the huge difference in the structural 
length scales of the two species. 



are of the order of the colloid size, not the much smaller star size. Altogether, these 
features show the tendency of colloids as well as stars to seek spatial proximity of their 
own species while sort of avoiding the other one and we may expect macroscopic regions 
rich in the one and poor in the other species to be formed provided the partial densities, 
in particular of the colloids, are sufficiently high. 

Fig. [7] illustrates the typical shape of the partial structure factors Sij(k). Here, 
we chose the parameters as follows: the PE-star functionality is f — 18, we set the 
size ratio to q — 0.2, fixed the density of the stars as p s a^ = 0.12, and considered 
several values of the colloidal density p c a^ . When comparing the three main plots of the 
figure, the first finding is that the locations of the different Lifshitz lines in density space 
strongly vary. These lines mark the respective structure factors' cross-over between a 
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Figure 8. Comparison of the cross structure factors S cs (k) for / = 30 and the two 
different size ratios investigated, q — 0.2 and q — 0.3. The PE-star partial densities 
were chosen to be almost the same in both cases, i.e., p s of = 0.27 and p s a^ = 0.24, 
respectively. The corresponding colloid densities are p c a^ = 0.29 (main plot) and 
p c a^ = 0.05 (inset). Upon varying the size ratio, the peak positions shift and the 
(k = 0)-values of the partial structure factors shown change significantly while there 
is no remarkable effect on the different peaks' height. 



regime where they display a local minimum in the long wavelength limit k — > and 
a region where the behaviour changes to developing a local maximum for the same k- 
values. While for the given amount of stars in the system the star-star Lifshitz line is 
obviously immediately crossed for practically arbitrary low colloid concentrations [fig. 
[7(b)], we need an noticeably increased partial colloid density lying in the range of about 
p c o~c ~ 0.25 . . . 0.5 for the colloid-star structure factor to experience such cross-over 
[fig- 0(c)]. In case of the colloid-colloid structure factor, the corresponding values of the 
colloid density are even higher, about p c o~l ~ 0.7 for the parameters used here [fig. Eta)]. 
Another indication of the demixing transition we are searching for within the scope of 
this paper and that is expected to occur upon adding more and more colloids and stars 
to the binary mixture is the tendency of all partial structure factors to diverge in the 
aforementioned long wavelength limit, i.e., S cc {k — ► 0) — ► +oo, S ss {k — > 0) — > +oo 
and S cs (k — ■> 0) —>■ — oo, thus demonstrating that we approach the spinodal line. The 
inset in fig. [7J(b) was included in order to again demonstrate the huge difference in 
the structural length scales of the two species present in the mixture. The pre-peak 
in the cross structure factor S cs (k) is without any direct physical interpretation, while 
pre-peaks in the intra-species structure factors would evince micro-phase separation 
[UJ SHI SSI [50] [51]. Since the latter peaks are completely absent in our case, we may 
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Figure 9. Demixing binodals calculated according to the procedures introduced 
in the main text for (a) q = 0.2 and (b) q = 0.3, and different values of the PE- 
stars' functionality /. In order to illustrate the coexisting colloid-poor and colloid-rich 
phases, we additionally show several tie lines. In this connection, please note that we 
in fact used much more such point pairs in order to obtain the binodal lines and not 
only the shown ones. Based on the full sets of coexisting fluid phases we computed, 
we made rough estimates for the positions of the respective critical points in the p c -p c 
plane, represented by the filled triangles. 



once more conclude that the system is expected to macro-phase separate instead of 
forming clusters. 

Finally, fig. [8] depicts the g-dependence of the cross structure factors for / = 30 
and two different values of the colloid density p c a^ (main plot and inset). For both 
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size ratios investigated, the star densities p s a^ are chosen to be almost the sam fl As 
obvious from the plots, a change in q only affects the peak positions and the depth of 
the local minimum for & — » 0, but there is no significant influence on the peak heights 
of the functions. This is in agreement with the findings for the g-dependence of the 
effective colloid-colloid interactions, see fig. El and essentially means that the size ratio 
q is crucial for determining the typical structural length scales, but hardly for how 
pronounced this structure is. 

4-3. Fluid- fluid phase equilibria 

After having found plenty of evidence in our hitherto analysis for a mixing-demixing 
transition taking place for certain ranges of partial densities Pi<jf, we finally come to a 
more quantitative description based on the corresponding binodals obtained as explained 
above. In fig. [9]we show the obtained demixing binodals for size ratios q = 0.2 [fig. [9](a)] 
and q = 0.3 [fig. Mjo)} and for different PE-star functionalities /, as denoted in the 
legend boxes. In addition, we connected some of the point pairs used to compute the 
binodals and representing coexisting colloid-rich and colloid-poor phases by tie lines. 
Concerning the mutual positions of the binodals in the density plane, it can be seen 
that they shift towards higher PE-star concentrations upon increasing the size ratio 
q and/or decreasing the PE-stars' functionality /. This characteristic behaviour is in 
agreement with previous studies of binary mixtures of colloids and neutral polymer stars 
[33] . The filled triangles in fig.[9]denote rough estimates for the respective critical points' 
positions determined graphically by taking the tie lines into account. The critical points 
move towards slightly lower colloid densities when lowering the PE-stars' arm number, 
whereas there is no significant effect of altering the size ratio. 

The star densities p s a^ that bring about a demixing instability are typically higher 
for the case q = 0.3 than for the case q = 0.2. This looks counterintuitive at first sight, 
since one expects that larger PE-stars will destabilise the mixture earlier. In order to 
put the numbers in their appropriate context, it is useful to employ the picture of the 
effective colloid-colloid potential, which includes a star-induced attraction. Here, the 
range and depth of this attraction steer the occurrence of the demixing binodal, which 
is equivalent to a separation between a colloidal fluid and a colloidal gas. The natural 
length scale in this picture is the colloid diameter a c ; concomitantly, the physically 
relevant density in making comparisons between the q = 0.2 and the q = 0.3 cases 
should be scaled with the colloid size: p s a^ = q~ 3 p s a 3 . It can be easily seen that the 
additional prefactor g -3 renders the rescaled star densities for q = 0.3 indeed lower than 
the ones for q = 0.2, in agreement with the intuitive expectations. 

The volume terms for the integrated out counterions [551 EEJ EI] do not affect the 
phase boundaries, since, under the assumption of full absorbing in the stars' interiors, 

* They are not exactly the same since such results are not systematically available due to the fact 
that we originally solved the OZ relation together with the RY closure for points in the density plane 
where the star density takes 'smooth' values when scaled with respect to the colloidal diameter a c , not 
their own diameter er s . 
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they are simply proportional to the number N s of the latter [TT] and thus they cause 
a trivial shift of the stars' chemical potential, without affecting the solution's osmotic 
pressure [58]. Finally, we mention that we did not consider the competition between 
the demixing binodals and the crystallization of the colloids. The investigation of 
the system's solid states lies beyond the scope of this work. The trends found for 
the /- and g-dependences are comparable to the colloid-star polymer case mentioned 
above. Although the underlying pair potentials are different to a certain degree, a closer 
inspection to the full phase diagrams in Ref. [33J can give hints regarding the stability 
of the binodals against preemption by the freezing lines. Provided the positions of the 
freezing lines are not too different here, it seems to be reasonable to assume based on 
such a comparison that our demixing lines will survive at least for the larger size ratio 
between stars and colloids. Nevertheless, the existence of a demixing binodal, even in 
the case that the latter is preempted by crystallization, has important consequences for 
the time scales involved in the dynamics of crystallization J59J [60] . 

5. Summary and conclusions 

We have put forward a coarse-grained description of mixtures between neutral, spherical, 
hard colloids and multiarm polyelectrolyte stars of size smaller than the colloidal 
particles. Effective interactions between the constituent particles have been employed 
throughout, allowing for a mesoscopic description that leads to valuable information on 
the structure and thermodynamics of the two-component mixture. The cross interaction, 
which has been derived in this work, is sufficiently repulsive to bring about regions of 
instability in the phase diagram and leading thereby to macroscopic, demixing phase 
behaviour. This, in turn, can be rationalised by means of the depletion potentials 
between the colloids, which are induced by the stars, and feature attractive tails akin 
to those encountered in usual colloid-polymer mixtures. 

The form of the cross interaction plays a crucial role in determining stability 
and can, by suitable tuning, completely change the behaviour of the mixture from 
macroscopic phase separation to microphase structuring with a finite wavelength. In 
this respect, a very promising direction of investigation is to allow for the colloids to 
carry a charge opposite to that of the arms of the polyelectrolyte stars. Preliminary 
results already indicate a rich variety of resulting complexation morphologies between 
the two constituents [61]. A detailed investigation of the complexation characteristics 
and the morphologies of the ensuing macroscopic phases is the subject of ongoing work. 
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